using DataFrames, Plots, Statistics, JSON, Distributed
plotlyjs();
The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.
The problem is simulated on 14 2.20 GHz processors running in parallel
const numcores = 14
if nprocs() < numcores
addprocs(numcores - nprocs())
end
@show nprocs();
@everywhere include(joinpath(dirname(pwd()),"jgrc/TuLiPa/src/TuLiPa.jl"))
@everywhere include("JulES/src/JulES.jl");
nprocs() = 14
datayear = 2025
scenarioyearstart = 1981
scenarioyearstop = 2011
horizonstart = 1981
numscen = 7
tnormal = FixedDataTwoTime(getisoyearstart(datayear),getisoyearstart(horizonstart))
phaseinoffsetdays = 2
phaseinoffset = Millisecond(Day(phaseinoffsetdays)) # phase in straight away from second stage scenarios
phaseindelta = Millisecond(Week(5)) # Phase in the second stage scenario over 5 weeks
phaseinsteps = 5 # Phase in second stage scenario in 5 steps
tphasein = PhaseinTwoTime(getisoyearstart(datayear),getisoyearstart(horizonstart), getisoyearstart(horizonstart), phaseinoffset, phaseindelta, phaseinsteps);
This is a simplified version of the dataset NVE uses for its long-term power market analyses (for example https://www.nve.no/energi/analyser-og-statistikk/langsiktig-kraftmarkedsanalyse/). The dataset consist of:
We cannot publish the full dataset, but many of the profiles for wind, solar, hydro and demand can be downloaded from https://www.nve.no/energi/analyser-og-statistikk/vaerdatasett-for-kraftsystemmodellene/.
sti_themadata = "data_fra_thema"
themastructure = json_to_elements(sti_themadata, "dataset_thema.json")
themaseries = json_to_elements(sti_themadata, "tidsserier_thema.json")
elements = vcat(themaseries,themastructure)
addscenariotimeperiod!(elements, scenarioyearstart, scenarioyearstop);
# Set horizons for price prognosis models
# Long
longhorizonend = getisoyearstart(horizonstart + 5)
longhydroperiodduration = Millisecond(Week(6))
longpowerparts = 6
longrhsdata = StaticRHSAHData("Power", datayear, scenarioyearstart, scenarioyearstop)
longmethod = KMeansAHMethod()
longclusters = 4
longunitduration = Millisecond(Hour(6))
# Medium
medweeklength = 54
medhorizonend = getisoyearstart(horizonstart) + Week(medweeklength)
medhydroperiodduration = Millisecond(Day(7)); @assert medweeklength % Int(longhydroperiodduration.value/3600000/168) == 0
medpowerparts = 7
medrhsdata = DynamicRHSAHData("Power")
medmethod = KMeansAHMethod()
medclusters = 4
medunitduration = Millisecond(Hour(4))
# Short
shortweeklength = 1
shorthorizonend = getisoyearstart(horizonstart) + Week(shortweeklength)
shorthydroperiodduration = Millisecond(Day(1)); @assert medweeklength % shortweeklength == 0
shortpowerparts = 12
# Preallocate results, distributed on cores
longprobs = distribute([HiGHS_Prob() for i in 1:numscen])
medprobs = distribute([HiGHS_Prob() for i in 1:numscen], longprobs)
shortprobs = distribute([HiGHS_Prob() for i in 1:numscen], longprobs)
medprices = distribute([Dict() for i in 1:numscen], longprobs)
shortprices = distribute([Dict() for i in 1:numscen], longprobs)
medendvaluesobjs = distribute([EndValues() for i in 1:numscen], longprobs)
nonstoragestates = distribute([Dict{StateVariableInfo, Float64}() for i in 1:numscen], longprobs)
# Organise inputs and outputs
probs = (longprobs, medprobs, shortprobs)
longinput = (longhorizonend, longhydroperiodduration, longrhsdata, longmethod, longclusters, longunitduration, "long")
medinput = (medhorizonend, medhydroperiodduration, medrhsdata, medmethod, medclusters, medunitduration, "long")
shortinput = (shorthorizonend, shorthydroperiodduration, shortpowerparts, "short")
allinput = (numcores, elements, horizonstart, tnormal, tphasein, phaseinoffsetdays)
output = (medprices, shortprices, medendvaluesobjs, nonstoragestates)
# Initialize price prognosis models and run for first time step. Run scenarios in parallell
@time pl_prognosis_init!(probs, allinput, longinput, medinput, shortinput, output)
43.805738 seconds (3.75 M allocations: 184.031 MiB, 0.20% gc time, 2.74% compilation time)
Task (done) @0x0000000098f73e10
scenario = 1
index = collect(medprices[scenario]["steprange"])
prices = medprices[scenario]["matrix"]
labels = [name for name in medprices[scenario]["names"]]
p = plot(index,prices,label=reshape(labels, 1, length(labels)),legend=:outertopright)
for scenario in 2:numscen
prices = medprices[scenario]["matrix"]
labels = [name for name in medprices[scenario]["names"]]
plot!(p,index,prices,label=reshape(labels, 1, length(labels)),legend=:outertopright)
end
display(p)